1   /*
2    * Copyright (C) 2011 The Guava Authors
3    *
4    * Licensed under the Apache License, Version 2.0 (the "License");
5    * you may not use this file except in compliance with the License.
6    * You may obtain a copy of the License at
7    *
8    * http://www.apache.org/licenses/LICENSE-2.0
9    *
10   * Unless required by applicable law or agreed to in writing, software
11   * distributed under the License is distributed on an "AS IS" BASIS,
12   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13   * See the License for the specific language governing permissions and
14   * limitations under the License.
15   */
16  
17  package com.google.common.math;
18  
19  import static com.google.common.base.Preconditions.checkArgument;
20  import static com.google.common.math.DoubleUtils.IMPLICIT_BIT;
21  import static com.google.common.math.DoubleUtils.SIGNIFICAND_BITS;
22  import static com.google.common.math.DoubleUtils.getSignificand;
23  import static com.google.common.math.DoubleUtils.isFinite;
24  import static com.google.common.math.DoubleUtils.isNormal;
25  import static com.google.common.math.DoubleUtils.scaleNormalize;
26  import static com.google.common.math.MathPreconditions.checkInRange;
27  import static com.google.common.math.MathPreconditions.checkNonNegative;
28  import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
29  import static java.lang.Math.abs;
30  import static java.lang.Math.copySign;
31  import static java.lang.Math.getExponent;
32  import static java.lang.Math.log;
33  import static java.lang.Math.rint;
34  
35  import com.google.common.annotations.GwtCompatible;
36  import com.google.common.annotations.GwtIncompatible;
37  import com.google.common.annotations.VisibleForTesting;
38  import com.google.common.primitives.Booleans;
39  
40  import java.math.BigInteger;
41  import java.math.RoundingMode;
42  import java.util.Iterator;
43  
44  /**
45   * A class for arithmetic on doubles that is not covered by {@link java.lang.Math}.
46   *
47   * @author Louis Wasserman
48   * @since 11.0
49   */
50  @GwtCompatible(emulated = true)
51  public final class DoubleMath {
52    /*
53     * This method returns a value y such that rounding y DOWN (towards zero) gives the same result
54     * as rounding x according to the specified mode.
55     */
56    @GwtIncompatible("#isMathematicalInteger, com.google.common.math.DoubleUtils")
57    static double roundIntermediate(double x, RoundingMode mode) {
58      if (!isFinite(x)) {
59        throw new ArithmeticException("input is infinite or NaN");
60      }
61      switch (mode) {
62        case UNNECESSARY:
63          checkRoundingUnnecessary(isMathematicalInteger(x));
64          return x;
65  
66        case FLOOR:
67          if (x >= 0.0 || isMathematicalInteger(x)) {
68            return x;
69          } else {
70            return x - 1.0;
71          }
72  
73        case CEILING:
74          if (x <= 0.0 || isMathematicalInteger(x)) {
75            return x;
76          } else {
77            return x + 1.0;
78          }
79  
80        case DOWN:
81          return x;
82  
83        case UP:
84          if (isMathematicalInteger(x)) {
85            return x;
86          } else {
87            return x + Math.copySign(1.0, x);
88          }
89  
90        case HALF_EVEN:
91          return rint(x);
92  
93        case HALF_UP: {
94          double z = rint(x);
95          if (abs(x - z) == 0.5) {
96            return x + copySign(0.5, x);
97          } else {
98            return z;
99          }
100       }
101 
102       case HALF_DOWN: {
103         double z = rint(x);
104         if (abs(x - z) == 0.5) {
105           return x;
106         } else {
107           return z;
108         }
109       }
110 
111       default:
112         throw new AssertionError();
113     }
114   }
115 
116   /**
117    * Returns the {@code int} value that is equal to {@code x} rounded with the specified rounding
118    * mode, if possible.
119    *
120    * @throws ArithmeticException if
121    *         <ul>
122    *         <li>{@code x} is infinite or NaN
123    *         <li>{@code x}, after being rounded to a mathematical integer using the specified
124    *         rounding mode, is either less than {@code Integer.MIN_VALUE} or greater than {@code
125    *         Integer.MAX_VALUE}
126    *         <li>{@code x} is not a mathematical integer and {@code mode} is
127    *         {@link RoundingMode#UNNECESSARY}
128    *         </ul>
129    */
130   @GwtIncompatible("#roundIntermediate")
131   public static int roundToInt(double x, RoundingMode mode) {
132     double z = roundIntermediate(x, mode);
133     checkInRange(z > MIN_INT_AS_DOUBLE - 1.0 & z < MAX_INT_AS_DOUBLE + 1.0);
134     return (int) z;
135   }
136 
137   private static final double MIN_INT_AS_DOUBLE = -0x1p31;
138   private static final double MAX_INT_AS_DOUBLE = 0x1p31 - 1.0;
139 
140   /**
141    * Returns the {@code long} value that is equal to {@code x} rounded with the specified rounding
142    * mode, if possible.
143    *
144    * @throws ArithmeticException if
145    *         <ul>
146    *         <li>{@code x} is infinite or NaN
147    *         <li>{@code x}, after being rounded to a mathematical integer using the specified
148    *         rounding mode, is either less than {@code Long.MIN_VALUE} or greater than {@code
149    *         Long.MAX_VALUE}
150    *         <li>{@code x} is not a mathematical integer and {@code mode} is
151    *         {@link RoundingMode#UNNECESSARY}
152    *         </ul>
153    */
154   @GwtIncompatible("#roundIntermediate")
155   public static long roundToLong(double x, RoundingMode mode) {
156     double z = roundIntermediate(x, mode);
157     checkInRange(MIN_LONG_AS_DOUBLE - z < 1.0 & z < MAX_LONG_AS_DOUBLE_PLUS_ONE);
158     return (long) z;
159   }
160 
161   private static final double MIN_LONG_AS_DOUBLE = -0x1p63;
162   /*
163    * We cannot store Long.MAX_VALUE as a double without losing precision.  Instead, we store
164    * Long.MAX_VALUE + 1 == -Long.MIN_VALUE, and then offset all comparisons by 1.
165    */
166   private static final double MAX_LONG_AS_DOUBLE_PLUS_ONE = 0x1p63;
167 
168   /**
169    * Returns the {@code BigInteger} value that is equal to {@code x} rounded with the specified
170    * rounding mode, if possible.
171    *
172    * @throws ArithmeticException if
173    *         <ul>
174    *         <li>{@code x} is infinite or NaN
175    *         <li>{@code x} is not a mathematical integer and {@code mode} is
176    *         {@link RoundingMode#UNNECESSARY}
177    *         </ul>
178    */
179   @GwtIncompatible("#roundIntermediate, java.lang.Math.getExponent, "
180       + "com.google.common.math.DoubleUtils")
181   public static BigInteger roundToBigInteger(double x, RoundingMode mode) {
182     x = roundIntermediate(x, mode);
183     if (MIN_LONG_AS_DOUBLE - x < 1.0 & x < MAX_LONG_AS_DOUBLE_PLUS_ONE) {
184       return BigInteger.valueOf((long) x);
185     }
186     int exponent = getExponent(x);
187     long significand = getSignificand(x);
188     BigInteger result = BigInteger.valueOf(significand).shiftLeft(exponent - SIGNIFICAND_BITS);
189     return (x < 0) ? result.negate() : result;
190   }
191 
192   /**
193    * Returns {@code true} if {@code x} is exactly equal to {@code 2^k} for some finite integer
194    * {@code k}.
195    */
196   @GwtIncompatible("com.google.common.math.DoubleUtils")
197   public static boolean isPowerOfTwo(double x) {
198     return x > 0.0 && isFinite(x) && LongMath.isPowerOfTwo(getSignificand(x));
199   }
200 
201   /**
202    * Returns the base 2 logarithm of a double value.
203    *
204    * <p>Special cases:
205    * <ul>
206    * <li>If {@code x} is NaN or less than zero, the result is NaN.
207    * <li>If {@code x} is positive infinity, the result is positive infinity.
208    * <li>If {@code x} is positive or negative zero, the result is negative infinity.
209    * </ul>
210    *
211    * <p>The computed result is within 1 ulp of the exact result.
212    *
213    * <p>If the result of this method will be immediately rounded to an {@code int},
214    * {@link #log2(double, RoundingMode)} is faster.
215    */
216   public static double log2(double x) {
217     return log(x) / LN_2; // surprisingly within 1 ulp according to tests
218   }
219 
220   private static final double LN_2 = log(2);
221 
222   /**
223    * Returns the base 2 logarithm of a double value, rounded with the specified rounding mode to an
224    * {@code int}.
225    *
226    * <p>Regardless of the rounding mode, this is faster than {@code (int) log2(x)}.
227    *
228    * @throws IllegalArgumentException if {@code x <= 0.0}, {@code x} is NaN, or {@code x} is
229    *         infinite
230    */
231   @GwtIncompatible("java.lang.Math.getExponent, com.google.common.math.DoubleUtils")
232   @SuppressWarnings("fallthrough")
233   public static int log2(double x, RoundingMode mode) {
234     checkArgument(x > 0.0 && isFinite(x), "x must be positive and finite");
235     int exponent = getExponent(x);
236     if (!isNormal(x)) {
237       return log2(x * IMPLICIT_BIT, mode) - SIGNIFICAND_BITS;
238       // Do the calculation on a normal value.
239     }
240     // x is positive, finite, and normal
241     boolean increment;
242     switch (mode) {
243       case UNNECESSARY:
244         checkRoundingUnnecessary(isPowerOfTwo(x));
245         // fall through
246       case FLOOR:
247         increment = false;
248         break;
249       case CEILING:
250         increment = !isPowerOfTwo(x);
251         break;
252       case DOWN:
253         increment = exponent < 0 & !isPowerOfTwo(x);
254         break;
255       case UP:
256         increment = exponent >= 0 & !isPowerOfTwo(x);
257         break;
258       case HALF_DOWN:
259       case HALF_EVEN:
260       case HALF_UP:
261         double xScaled = scaleNormalize(x);
262         // sqrt(2) is irrational, and the spec is relative to the "exact numerical result,"
263         // so log2(x) is never exactly exponent + 0.5.
264         increment = (xScaled * xScaled) > 2.0;
265         break;
266       default:
267         throw new AssertionError();
268     }
269     return increment ? exponent + 1 : exponent;
270   }
271 
272   /**
273    * Returns {@code true} if {@code x} represents a mathematical integer.
274    *
275    * <p>This is equivalent to, but not necessarily implemented as, the expression {@code
276    * !Double.isNaN(x) && !Double.isInfinite(x) && x == Math.rint(x)}.
277    */
278   @GwtIncompatible("java.lang.Math.getExponent, com.google.common.math.DoubleUtils")
279   public static boolean isMathematicalInteger(double x) {
280     return isFinite(x)
281         && (x == 0.0 ||
282             SIGNIFICAND_BITS - Long.numberOfTrailingZeros(getSignificand(x)) <= getExponent(x));
283   }
284 
285   /**
286    * Returns {@code n!}, that is, the product of the first {@code n} positive
287    * integers, {@code 1} if {@code n == 0}, or {@code n!}, or
288    * {@link Double#POSITIVE_INFINITY} if {@code n! > Double.MAX_VALUE}.
289    *
290    * <p>The result is within 1 ulp of the true value.
291    *
292    * @throws IllegalArgumentException if {@code n < 0}
293    */
294   public static double factorial(int n) {
295     checkNonNegative("n", n);
296     if (n > MAX_FACTORIAL) {
297       return Double.POSITIVE_INFINITY;
298     } else {
299       // Multiplying the last (n & 0xf) values into their own accumulator gives a more accurate
300       // result than multiplying by everySixteenthFactorial[n >> 4] directly.
301       double accum = 1.0;
302       for (int i = 1 + (n & ~0xf); i <= n; i++) {
303         accum *= i;
304       }
305       return accum * everySixteenthFactorial[n >> 4];
306     }
307   }
308 
309   @VisibleForTesting
310   static final int MAX_FACTORIAL = 170;
311 
312   @VisibleForTesting
313   static final double[] everySixteenthFactorial = {
314       0x1.0p0,
315       0x1.30777758p44,
316       0x1.956ad0aae33a4p117,
317       0x1.ee69a78d72cb6p202,
318       0x1.fe478ee34844ap295,
319       0x1.c619094edabffp394,
320       0x1.3638dd7bd6347p498,
321       0x1.7cac197cfe503p605,
322       0x1.1e5dfc140e1e5p716,
323       0x1.8ce85fadb707ep829,
324       0x1.95d5f3d928edep945};
325 
326   /**
327    * Returns {@code true} if {@code a} and {@code b} are within {@code tolerance} of each other.
328    *
329    * <p>Technically speaking, this is equivalent to
330    * {@code Math.abs(a - b) <= tolerance || Double.valueOf(a).equals(Double.valueOf(b))}.
331    *
332    * <p>Notable special cases include:
333    * <ul>
334    * <li>All NaNs are fuzzily equal.
335    * <li>If {@code a == b}, then {@code a} and {@code b} are always fuzzily equal.
336    * <li>Positive and negative zero are always fuzzily equal.
337    * <li>If {@code tolerance} is zero, and neither {@code a} nor {@code b} is NaN, then
338    * {@code a} and {@code b} are fuzzily equal if and only if {@code a == b}.
339    * <li>With {@link Double#POSITIVE_INFINITY} tolerance, all non-NaN values are fuzzily equal.
340    * <li>With finite tolerance, {@code Double.POSITIVE_INFINITY} and {@code
341    * Double.NEGATIVE_INFINITY} are fuzzily equal only to themselves.
342    * </li>
343    *
344    * <p>This is reflexive and symmetric, but <em>not</em> transitive, so it is <em>not</em> an
345    * equivalence relation and <em>not</em> suitable for use in {@link Object#equals}
346    * implementations.
347    *
348    * @throws IllegalArgumentException if {@code tolerance} is {@code < 0} or NaN
349    * @since 13.0
350    */
351   public static boolean fuzzyEquals(double a, double b, double tolerance) {
352     MathPreconditions.checkNonNegative("tolerance", tolerance);
353     return
354           Math.copySign(a - b, 1.0) <= tolerance
355            // copySign(x, 1.0) is a branch-free version of abs(x), but with different NaN semantics
356           || (a == b) // needed to ensure that infinities equal themselves
357           || (Double.isNaN(a) && Double.isNaN(b));
358   }
359 
360   /**
361    * Compares {@code a} and {@code b} "fuzzily," with a tolerance for nearly-equal values.
362    *
363    * <p>This method is equivalent to
364    * {@code fuzzyEquals(a, b, tolerance) ? 0 : Double.compare(a, b)}. In particular, like
365    * {@link Double#compare(double, double)}, it treats all NaN values as equal and greater than all
366    * other values (including {@link Double#POSITIVE_INFINITY}).
367    *
368    * <p>This is <em>not</em> a total ordering and is <em>not</em> suitable for use in
369    * {@link Comparable#compareTo} implementations.  In particular, it is not transitive.
370    *
371    * @throws IllegalArgumentException if {@code tolerance} is {@code < 0} or NaN
372    * @since 13.0
373    */
374   public static int fuzzyCompare(double a, double b, double tolerance) {
375     if (fuzzyEquals(a, b, tolerance)) {
376       return 0;
377     } else if (a < b) {
378       return -1;
379     } else if (a > b) {
380       return 1;
381     } else {
382       return Booleans.compare(Double.isNaN(a), Double.isNaN(b));
383     }
384   }
385 
386   @GwtIncompatible("com.google.common.math.DoubleUtils")
387   private static final class MeanAccumulator {
388 
389     private long count = 0;
390     private double mean = 0.0;
391 
392     void add(double value) {
393       checkArgument(isFinite(value));
394       ++count;
395       // Art of Computer Programming vol. 2, Knuth, 4.2.2, (15)
396       mean += (value - mean) / count;
397     }
398 
399     double mean() {
400       checkArgument(count > 0, "Cannot take mean of 0 values");
401       return mean;
402     }
403   }
404 
405   /**
406    * Returns the arithmetic mean of the values. There must be at least one value, and they must all
407    * be finite.
408    */
409   @GwtIncompatible("MeanAccumulator")
410   public static double mean(double... values) {
411     MeanAccumulator accumulator = new MeanAccumulator();
412     for (double value : values) {
413       accumulator.add(value);
414     }
415     return accumulator.mean();
416   }
417 
418   /**
419    * Returns the arithmetic mean of the values. There must be at least one value. The values will
420    * be converted to doubles, which does not cause any loss of precision for ints.
421    */
422   @GwtIncompatible("MeanAccumulator")
423   public static double mean(int... values) {
424     MeanAccumulator accumulator = new MeanAccumulator();
425     for (int value : values) {
426       accumulator.add(value);
427     }
428     return accumulator.mean();
429   }
430 
431   /**
432    * Returns the arithmetic mean of the values. There must be at least one value. The values will
433    * be converted to doubles, which causes loss of precision for longs of magnitude over 2^53
434    * (slightly over 9e15).
435    */
436   @GwtIncompatible("MeanAccumulator")
437   public static double mean(long... values) {
438     MeanAccumulator accumulator = new MeanAccumulator();
439     for (long value : values) {
440       accumulator.add(value);
441     }
442     return accumulator.mean();
443   }
444 
445   /**
446    * Returns the arithmetic mean of the values. There must be at least one value, and they must all
447    * be finite. The values will be converted to doubles, which may cause loss of precision for some
448    * numeric types.
449    */
450   @GwtIncompatible("MeanAccumulator")
451   public static double mean(Iterable<? extends Number> values) {
452     MeanAccumulator accumulator = new MeanAccumulator();
453     for (Number value : values) {
454       accumulator.add(value.doubleValue());
455     }
456     return accumulator.mean();
457   }
458 
459   /**
460    * Returns the arithmetic mean of the values. There must be at least one value, and they must all
461    * be finite. The values will be converted to doubles, which may cause loss of precision for some
462    * numeric types.
463    */
464   @GwtIncompatible("MeanAccumulator")
465   public static double mean(Iterator<? extends Number> values) {
466     MeanAccumulator accumulator = new MeanAccumulator();
467     while (values.hasNext()) {
468       accumulator.add(values.next().doubleValue());
469     }
470     return accumulator.mean();
471   }
472 
473   private DoubleMath() {}
474 }